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Cyclic dominance and biodiversity in well-mixed populations 
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Coevolutionary dynamics is investigated in chemical catalysis, biological evolution, social and economic sys- 
tems. The dynamics of these systems can be analyzed within the unifying framework of evolutionary game 
theory. In this Letter, we show that even in well-mixed finite populations, where the dynamics is inherently 
stochastic, biodiversity is possible with three cyclic dominant strategies. We show how the interplay of evolu- 
tionary dynamics, discreteness of the population, and the nature of the interactions influences the coexistence of 
strategies. We calculate a critical population size above which coexistence is likely. 

PACS numbers: 87.23.-n, 89.65.-s 02.50.Ey 



Coevolution with cyclic dominance can lead to oscillatory, 
chaotic and stochastic dynamics. For example, such cycles 
can be found in biology lilSsflSlsLZl] or in social and eco- 
nomic systems Jl, HI] . The simplest example for such a cyclic 
dominance is the well known children's game rock-paper- 
scissors, where rock crushes scissors, scissors cuts paper and 
paper wraps rock. Different biological realizations of this sys- 
tem have been observed. For example, a cyclic dominance 
of three male strategies has been reported in lizards 01 0]: 
Orange-throated males establish large territories holding sev- 
eral females. These populations are invaded by males with 
yellow-striped throats, which do not contribute to the defense 
of the territory but sneak on the females. Such a population of 
yellow-striped males can be invaded by blue-throated males, 
which defend territories large enough to hold one female 
which they defend against sneakers. Once yellow-striped 
sneakers are rare, it is advantageous to defend a large territory 
with several females and the cycle starts with orange-throated 
males again. Another example is the competition between dif- 
ferent strains of E.coli. Kerr et aJ. [3, 4] observed that cyclic 
dominance leads to biodiversity in spatial systems, whereas 
two strategies go extinct in mixed systems after a short time. 
Such systems have been analyzed by evolutionary game the- 
ory in great mathematical detail 11 1 QL [HI El [II Hill a. 
Based on the replicator dynamics describing the dynamics in 
an infinite population, general conclusions on the nature of 
the payoff matrix can be made from the observation of fast 
extinction. How does this picture change if we tackle the 
more realistic case of stochastic dynamics in a finite popu- 
lation 11171 Il8l Il9ll ? Dependin g on popula tion size and the 
underlying microscopic process l20l |2 ll l22l l23ll . the resulting 
dynamics can be very different from the the replicator equa- 
tion results. 



Definition of the model. - We first concentrate on the 
simplest example of cyclic rock-paper-scissors dynamics, in 
which all three strategies are equivalent [24]. Thus, the game 
is defined by three payoffs: (i) the payoff against a dominated 
strategy (set to 1). (ii) the payoff against a dominant strategy 
— s, which we assume to be negative, (iii) the payoff for a tie 



(set to 0). Thus, we obtain the payoff matrix 
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Only for the standard choice s = 1, we have a zero-sum game. 
An intuitive understanding how s influences the game can be 
obtained from discussing two relevant cases. Large values of s 
will make it successful to avoid losing, best done with staying 
with the majority. In this case, a mixed equilibrium is unstable 
and ultimately, only one strategy will survive. For s f=s 0, it 
is more important to win occasionally, such that the mixed 
equilibrium can become stable. 

Let us first recall the evolutionary dynamics in an infinite 
population described by the replicator dynamics [10]. In the 
replicator equation, the frequency (abundance) Xk of strategy 
k changes proportional to its payoff Hk, 



Xk = X k {TTk ~ <7r)) • 



(2) 



Here (ir) is the average payoff in the population. We denote 
the frequencies of R, P, S as x,y,z, respectively, with x + 
y + z — 1. With the standard assumption that payoffs are 
determined from interactions with a representative subset of 
the population, we find itr : = z — sy, up = x — sz, and tts — 
y — sx. The average payoff is given by (ir) = (1 — s)(xy + 
xz + yz) vanishing for the zero-sum game with s = 1. Eq. (0 
has an interior fixed point p = (| , i , | ) . The determinant of 
the payoff matrix d = 1 — s 3 determines the dynamics of the 
system [ 10]. For d > 0, which is the case for s < 1, 
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is a Lyapunov function with H < and the interior fixed point 
p is asymptotically stable. For d < (s > 1), p is unstable 
and the attractor of the system approaches a heteroclinic cycle 
at the boundary of the simplex S3. Finally, for the zero-sum 
game with s = 1 and d = the function (O is a constant of 
motion, and the system infinitely oscillates around p. This is 
not a purely theoretical exercise: Experiments indicate d > 
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for the Lizard system [2] and d < for the E.coli system 
0]. In theoretical approaches, often the restriction to s = 1 
is made, although this is a special non-generic case lfl5i[l6ll . 
It has been argued that in this case, limited mobility in a spa- 
tially extended system can promote biodiversity. However, for 
s < 1 biodiversity is likely even in a well-mixed population 
if the population is larger than N c . In this Letter, we give an 
analytical estimate for this critical population size, depending 
on s and the update mechanism. For s = 1 the average time 
to extinction scales linearly with the population size N 1)2311 - 
For s > 1, one expects faster extinction, as the stochastic drift 
and the instability of the interior fixed point act jointly. In the 
most interesting case of s < 1, both forces work in opposite 
directions: For N — > oo, the fixed point p is stable, but a finite 
N leads to a drift towards the simplex boundaries. 

Evolutionary processes. - To study dynamics in finite pop- 
ulations, we study microscopic stochastic processes for the 
dynamics, which lead to macroscopic equations of motion for 
large populations (|8l |20l |2l||. As the microscopic dynamics 
may depend on the system, the respective biological or be- 
havioral setup may require different interaction and competi- 
tion processes. To demonstrate the robustness of our results, 
we consider different birth-death processes: the frequency- 
dependent Moran process (MO) 17 1 ll &JlSffl, an d local two- 
particle interaction processes Il20l l2Ul27Ll28l 12911 . 

In finite populations with i i?-players, j P-players, and 
N — i — j S-players, the payoffs can be calculated from 
the equations for infinite N by replacing x — > i/(N — 1), 
V -> j/(N - 1), and z -> (N - i - j)/{N - 1). By di- 
viding by N — 1 and setting the payoff for ties zero, we for- 
mally exclude self-interactions. For the average payoff, we 



have (tt) = 4m R + U p + ^Mtt s 
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In the frequency-dependent Moran process, an individual 
reproduces proportional to its fitness. Then, the offspring re- 
places a randomly selected individual [30]. The transition 
probabilities of the possible six hopping events are given by 
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(4) 
(5) 



and T SP \T PS ',T PR ,T RP are obtained by cyclic permuta- 



tion of (R,P,S) and (i,j,N — i — j). Fitness is a convex 
combination of a constant background fitness (set to 1) and 
the payoff. The parameter w > controls the intensity of se- 
lection; random genetic drift is obtained for w — » 0. For better 
comparison with the processes below, we have introduced an 
additional factor 1/2. 

Selecting an individual proportional to fitness requires 
knowledge about every payoff in the population. In many 
cases, it is more realistic to assume that competition occurs 
locally between two individuals. One process of this type is 
the local update (LU) process |2(X 12111 . where one individ- 
ual b is selected randomly for reproduction, compares with 
annother randomly chosen individual a, and changes strategy 



with probability |(1 + w(n a — 7Tb)) 13111 . In general, the re 
productive fitness can depend in a nonlinear way on the pay- 
off difference between two competing agents, as in the Fermi 

process (FP) IHIH. 

We now unify the processes by means of a reproductive 
function, 

*Mo(o^a) = t,- ; — ( 6 ) 



2 1 — w + w{ti) 
$lm(6^o) = (1 + iw(7r„ - (tt»)/2 (7) 
$Lu(fc->a) = (1 + w(ir a - 7T 6 ))/2 (8) 
$ F p(fo->a) = [l + cxp(-ii;(7r a -7Tb)))]- 1 . (9) 
By (|7l) we introduce a linearized Moran (LM) process as first 
order approximation of the MO in the limit of weak selection, 
w — > 0. For w = 0, selection is neutral and the four processes 
are identical with random drift. Note that for the Moran 
processes, <i>(6 — * a) depends on the average payoff, whereas 
for the other processes two individual payoffs are involved. 
The transition probabilities become 



T ba = §{b-> a)N a N b / N 2 



(10) 



For all processes, $(6 — > a) considers a two-particle (birth- 
death) process where an individual with fitness 7r a compares 
with the average fitness (tt) or with annother individual Tib- 
Average drift. - For the replicator equation of the symmet- 
ric RPS dynamics, Eq. (01 defines a constant of motion. As we 
are interested in the finite-size corrections, we can use H as an 
observable for the distance to the interior fixed point. For the 
processes defined above, the transition probabilities allow to 
calculate the average change of the constant of motion within 
the simplex fl < j, j < N - 1) as ll32tl 



AT-lTV-i-l 

(AH)=^ E J2 E [ij(N-i-j)(T RS + T SR + T SP +T ps +T PR + T RP ) (11) 

»=1 3=1 

- (t - l)j(N -i-j+ l)T RS - (i + l)j(N -i-j- l)T SR - i(j + 1)(N — i — j — 1)T SP 

- i(J - 1)(N \)T PS - (i + -l)(N-i- j)T PR - (i - + 1)(N - i - j)T RP ^ 



These expressions give the exact average drift. For the two linear processes, we can approximate the average drift by replacing 
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the sums by integrals. Using x = i/N, y = j/N, and z — 1 — x — y, we find in the continuum limit 

^1 rl-x 



fl /■!— x 



y(x - z) (T RS - T SR ) + x(z - y) (T SP - T PS ) + z(y - x) (T 
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(12) 



From this expression, we can perform a comparison of (AH) 
for the different processes. The neutral case as well as the 
linear cases are analyzed below analytically, in Fig. 1 all pro- 
cesses are compared numerically. 
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FIG. 1: (Color online.) Drift reversal for the four processes; the 
Linearized Moran and Fermi process coincide within line thickness. 
The main panel shows the scaling with population size N for fixed 
w = 0.45 and s = 0.8. Full lines and the dotted line show the 
numerical solutions of Eq. d 1 1 b . The dashed lines are the analytical 
expansions Eqs. J 15 b and ( 1 17b . For small N, the continuum approx- 
imation to the integral leads to deviations. Inset: Average drift for 
varying payoff matrix with fixed w = 0.45 and N = 1000. The 
processes intersect in AH = l/(20iV 2 ) at s = 1 (d = 0) . For 
N — » oo, this reduces to the expectation based on the replicator 
equation, AH — in d — 0. The intersection with the horizontal 
line indicates the drift reversal in finite populations. 

For neutral evolution, w = 0, all terms of type T RS ~T SR 
vanish, and we have T RS = T SR = xz/2, T SP = T PS = 
yz/2, T PR = T RP = xy/2, i.e. for w = 0, 
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Now we consider the linear terms in w. For the linear local 
update, we find T RS - T SR = w xz(n s - n R ), thus 
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RS 



T SR = wxz((y - z) + s(y - x)). 



(14) 



Due to &Lu(b — > a) + $lu(o- —> b) = 1, the terms of order 
7V~ 2 in Eq. ( fT2l cancel. In first and second order of w, we 
have 

{AH)W - 2oW ~ (15) 



For the linearized Moran process, the contributions of the 
average payoff vanish for the drift term, and the calculation 
above reproduces with an additional factor 1/2, i.e. T RS — 
rpSR _ w xz f^ n s _ n Ry p or tne (jiff us i ori terrrli we nn( j 

t rs + t sr = aJ z(l + |!(7r sr + 7r JI -2(7r))) (16) 

and cyclic permutations. Thus, we get an additional contribu- 
tion to the diffusion. In summary, we have 



(AH) hM = 



1 1 - s ( 1 1 
20iV 2 ~ 42(W V2 ~ N 



(17) 




FIG. 2: (Color online.) Critical population size: Intensity of selection 
w versus critical population size iV c , where AH changes its sign 
for fixed s — 0.8. For N > N c the drift is towards the internal 
fixed point. The dashed lines refer to the Lizard's payoff matrix Q|] 
for MO and LU. Under strong selection w = 1 the LU reverses at 
N c — 20. For the restricted random matrix ensemble (green circles, 
see text) under LU N c coincides with the case of eq.fi] Each symbol 
corresponds to a single random payoff matrix, (inset: scatter plot of 
payoff entries in A 1 st , T 2 nd , ■ 3 rd column). 

Biodiversity threshold. - Both Eqs. <TT3T > and (T% can 
change their sign for s < 1 depending on N. For all four pro- 
cesses, Fig. [T]shows the average drift for different payoffs and 
population sizes. All processes intersect in +1/(2(W 2 ) for 
s = 1. In this case, we have a zero-sum game and the average 
drift is equal to neutral selection. The finite-size correction 
to the vanishing drift term of neutral selection arises from the 
difference between our microscopic processes and symplectic 
integrators ll25[ l26ll . In the — > oo limit we recover that the 



4 



stability of the interior fixed point is governed by the sign of 
the determinant of the payoff matrix, d = 1 — s 3 . For s < 1 
and N — > oc, trajectories spiral inwards. In finite populations, 
for s < 1 the stochastic motion can be reversed as shown in 
Fig- HI The critical population size of the biodiversity thresh- 
old is N c = -g—r for the LU and N c = 2 + -i^ for 

L w(l — s) L w(l — s) 

the LM. The biodiversity threshold N c for the other process is 
computed numerically in Fig. |2] For N > N c , the N — » oo 
behaviour is recovered, but for N < N c the interior fixed 
point becomes unstable and the average drift goes towards the 
boundaries. So far, we have worked with the specific payoff 
matrix Eq. ((T). Next, we show that this is a more general phe- 
nomenon. The drift reversal is preserved for general cyclic 
payoffs with N c depending on the payoff matrix [33]. As an 
example, we consider the payoff matrix of the Lizard system 
|Q]]. In this case, we find N c > 20, see Fig. [2] More general, 
we can also address random payoff matrices with cyclic dom- 
inance. Confining a random matrix ensemble (Fig. [U onto 
a suitable 2-dim submanifold [34], we show that N c (w) de- 
pends only on the location of the fixed point p, its normaliza- 
tion constant S and the determinant d, i.e. on 4 parameters. 
Hence the phenomenon is generic for arbitrary cyclic payoff 
matrices with d > 0. 



To conclude, cyclic coevolution in biological or social dy- 
namics highlights the importance to study finite population ef- 
fects and the underlying microscopic dynamics. Recently the 
influence of the finiteness of the population has been widely 
discussed for 2x2 games |Qj| [T^, HH H2]. But also in 
cyclic dynamics the finiteness of the population can modify 
the stability conditions derived from the replicator equation. 
In this Letter, we have shown that biodiversity threshold of the 
population size occurs in generic cyclic 3x3 games for a posi- 
tive determinant of the payoff matrix. Such a positive determi- 
nant has been found in the Lizard system |ji|2j]. Thus, stability 
of the coexistence fixed point can be obtained for sufficiently 
large populations, preserving biodiversity even in non-spatial 
systems or under strong mobility. In contrast, experiments in 
well-mixed populations of the E.coli system indicate a nega- 
tive determinant of the payoff matrix. 

Here we have demonstrated how nonzero-sum payoffs can 
change the stability even in well-mixed systems. Thus in bio- 
logical — and corresponding social or economic — systems, 
payoff matrix, population structure, population size and the 
microscopic update mechanism determine the fate of extinc- 
tion or coexistence. 



[1] B. Sinervo, CM. Lively, Nature 380, 240 (1996). 

[2] K.R. Zamudio and B. Sinervo, PNAS 97, 14427 (2000). 

[3] B. Kerr, M.A. Riley, M.W. Feldman, and B.J.M. Bohannan, Na- 
ture (London) 418, 171 (2002). 

[4] B.C. Kirkup and M.A. Riley, Nature 428, 412 (2004). 

[5] A. E. Magurran and M. A. Nowak, Proc. Roy. Soc. B 246, 31 
(1991). 

[6] Richard Dawkins, The Selfish Gene, Oxford UP (1976). 

[7] U. Bastolla, M. Lassig, S. C. Manrubia and A. Valleriani, J. 

Theor. Biol. 235, 521-530 (2005). 
[8] D. Helbing, Physica A 181, 29 (1992); 193, 241 (1993). 
[9] C. Hauert, S. De Monte, J. Hofbauer and K. Sigmund, Science 

296,1129 (2002). 
[10] J. Hofbauer, K. Sigmund, Evolutionary Games and Population 

Dynamics, Cambridge University Press (1998). 
[11] M. A. Nowak, Evolutionary Dynamics, Harvard University 

Press (2006). 

[12] G. Szabo and G. Fath, Phys. Rep. 446, 97 (2007). 
[13] G. Szabo and T. Czaran, Phys. Rev. E 64, 042902 (2001). 
[14] T. Czaran, R. F. Hoekstra, and L. Pagie, PNAS 99, 786 (2002). 
[15] T. Reichenbach, M. Mobilia, and E. Frey, Nature 448, 1046 
(2007). 

[16] T. Reichenbach, M. Mobilia, and E. Frey, Phys. Rev. Lett. 99, 
238105 (2007). 

[17] P. A. P. Moran, The Statistical Processes of Evolutionary The- 
ory, Clarendon, Oxford (1962). 

[18] M. A. Nowak, A. Sasaki, C. Taylor, and D. Fudenberg, Nature 
428, 646 (2004). 

[19] C. Taylor, D. Fudenberg, A. Sasaki, and M. A. Nowak, Bull. 
Math. Biol. 66, 1621 (2004). 



[20] A. Traulsen, J.C. Claussen, and C. Hauert, Phys. Rev. Lett. 95, 
238701 (2005). 

[21] A. Traulsen, J.C. Claussen, and C. Hauert, Phys. Rev. E 74, 
011901 (2006). 

[22] J.C. Claussen and A. Traulsen, Phys. Rev. E 71, 025 101 (R) 
(2005). 

[23] T. Reichenbach, M. Mobilia, and E. Frey, Phys. Rev. E 74, 
051907 (2006). 

[24] William H. Sandholm, Population Games and Evolutionary Dy- 
namics, MIT Press, in press. 

[25] J. Hofbauer, J. Math. Biol. 34, 675 (1996). 

[26] Y. Sato and J. P. Crutchfield, Phys. Rev. E 67, 015206(R) 
(2003). 

[27] L.E. Blume, Games Econ. Behav. 5, 387 (1993). 
[28] G. Szabo and C. Hauert, Phys. Rev. Lett. 89, 118101 (2002). 
[29] A. Traulsen, M. A. Nowak, and J. M. Pacheco, Phys. Rev. E 74, 
011909 (2006). 

[30] Reproduction proportional to fitness implies that 1 — w + wn 
must be positive. In our case, this leads to w < 1/2. 

[31] In Ref. 12011 , w is replaced by tu/A7r max . 

[32] Since we are only interested in large N, a prefactor ( N— 1) (N — 
2) /N 2 is omitted. 

[33] Using H := -P A with P = x Pj! y Py z Pz (where p x , p y and p z 
are the coordinates of the fixed point) and normalization con- 
stant E, see 1 10] for details. 

[34] With two payoff matrix entries equiprobable within a finite in- 
terval; remaining nondiagonal entries were computed such that 
the coordinates of the fixed point and E keep fixed values. 



